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Euclidean Distance Matrices 

Essential Theory, Algorithms and Applications 

Ivan Dokmanic, Reza Parhizkar, Juri Ranieri and Martin Vetterli 


Abstract —Euclidean distance matrices (EDM) are matrices of 
squared distances between points. The definition is deceivingly 
simple: thanks to their many useful properties they have found 
applications in psychometrics, crystallography, machine learn¬ 
ing, wireless sensor networks, acoustics, and more. Despite the 
usefulness of EDMs, they seem to be insufficiently known in the 
signal processing community. Our goal is to rectify this mishap 
in a concise tutorial. 

We review the fundamental properties of EDMs, such as 
rank or (non)definiteness. We show how various EDM properties 
can be used to design algorithms for completing and denoising 
distance data. Along the way, we demonstrate applications to 
microphone position calibration, ultrasound tomography, room 
reconstrnction from echoes and phase retrieval. By spelling out 
the essential algorithms, we hope to fast-track the readers in 
applying EDMs to their own problems. Matlab code for all the 
described algorithms, and to generate the figures in the paper, 
is available online. Finally, we suggest directions for further 
research. 

I. Introduction 

Imagine that you land at Geneva airport with the Swiss 
train schedule but no map. Perhaps surprisingly, this may be 
sufficient to reconstruct a rough (or not so rough) map of 
the Alpine country, even if the train times poorly translate to 
distances or some of the times are unknown. The way to do 
it is by using Euclidean distance matrices (EDM): for a quick 
illustration, take a look at the “Swiss Trains” box. 

An EDM is a matrix of squared Euclidean distances between 
points in a set.' We often work with distances because they 
are convenient to measure or estimate. In wireless sensor 
networks for example, the sensor nodes measure received 
signal strengths of the packets sent by other nodes, or time-of- 
arrival (TOA) of pulses emitted by their neighbors [1]. Both 
of these proxies allow for distance estimation between pairs 
of nodes, thus we can attempt to reconstruct the network 
topology. This is often termed self-localization [2], [3], [4]. 
The molecular conformation problem is another instance of a 
distance problem [5], and so is reconstructing a room’s geom¬ 
etry from echoes [6]. Less obviously, sparse phase retrieval [7] 
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Fig. 1. Two real-world applications of EDMs. Sensor network localization 
from estimated pairwise distances is illustrated on the left, with one distance 
missing because the con'esponding sensor nodes are too far apart to communi¬ 
cate. In the molecular conformation problem on the right, we aim to estimate 
the locations of the atoms in a molecule from their pairwise distances. Here, 
due to the inherent measurement uncertainty, we know the distances only up 
to an interval. 

can be converted to a distance problem, and addressed using 
EDMs. 

Sometimes the data are not metric, but we seek a metric 
representation, as it happens commonly in psychometrics [8]. 
As a matter of fact, the psychometrics community is at the root 
of the development of a number of tools related to EDMs, 
including multidimensional scaling (MDS)—the problem of 
finding the best point set representation of a given set of 
distances. More abstractly, people are concerned with EDMs 
for objects living in high-dimensional vector spaces, such as 
images [9]. 

EDMs are a useful description of the point sets, and a 
starting point for algorithm design. A typical task is to retrieve 
the original point configuration: it may initially come as a 
surprise that this requires no more than an eigenvalue decom¬ 
position (EVD) of a symmetric matrix.^ In fact, the majority 
of Euclidean distance problems require the reconstruction of 
the point set, but always with one or more of the following 
twists: 

1) Distances are noisy, 

2) Some distances are missing, 

^Because the EDMs are symmetric, we choose to use EVDs instead of 
singular value decompositions. That the EVD is much more efficient for 
symmetric matrices was suggested to us by one of the reviewers of the 
initial manuscript, who in turn received the advice from the numerical analyst 
Michael Saunders. 
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Swiss Trains (Swiss Map Reconstruction) 



Consider the following matrix of times in minutes it takes to 
travel by train between some Swiss cities: 



L 

G 

Z 

N 

B 

Lausanne 

f ° 

33 

128 

40 

66 \ 

Geneva 

33 

0 

158 

64 

101 

Zurich 

128 

158 

0 

88 

56 

Neuchatel 

40 

64 

88 

0 

34 

Bern 

\ 66 

101 

56 

34 

0 / 


The numbers were taken from the Swiss railways timetable. 
The matrix was then processed using the classical MDS algo¬ 
rithm (Algorithm 1), which is basically an EVD. The obtained 
city configuration was rotated and scaled to align with the 
actual map. Given all the uncertainties involved, the fit is 
remarkably good. Not all trains drive with the same speed; 
they have varying numbers of stops, railroads are not straight 
lines (lakes and mountains). This result may be regarded as 
anecdotal, but in a fun way it illustrates the power of the EDM 
toolbox. Classical MDS could be considered the simplest of the 
available tools, yet it yields usable results with erroneous data. 
On the other hand, it might be that Swiss trains are just so 
good. 


3) Distances are unlabeled. 

For two examples of applications requiring a solution of EDM 
problems with different complications, see Fig. 1. 

There are two fundamental problems associated with dis¬ 
tance geometry [10]: (i) given a matrix, determine whether it 
is an EDM, (ii) given a possibly incomplete set of distances, 
determine whether there exists a configuration of points in a 
given embedding dimension—dimension of the smallest affine 
space comprising the points—that generates the distances. 

A. Prior Art 

The study of point sets through pairwise distances, and so 
of EDMs, can be traced back to the works of Menger [11], 
Schoenberg [12], Blumenthal [13], and Young and House¬ 
holder [14]. 

An important class of EDM tools was initially developed 
for the purpose of data visualization. In 1952, Torgerson 
introduced the notion of MDS [8]. He used distances to quan¬ 
tify the dissimilarities between pairs of objects that are not 


necessarilly vectors in a metric space. Later in 1964, Kruskal 
suggested the notion of stress as a measure of goodness-of- 
ht for non-metric data [15], again representing experimental 
dissimilarities between objects. 

A number of analytical results on EDMs were developed 
by Gower [16], [17]. In his 1985 paper [17], he gave a 
complete characterization of the EDM rank. Optimization with 
EDMs requires good geometric intuitions about matrix spaces. 
In 1990, Glunt [18] and Hayden [19] with their co-authors 
provided insights into the structure of the convex cone of 
EDMs. An extensive treatise on EDMs with many original 
results and an elegant characterization of the EDM cone is 
given by Dattorro [20]. 

In early 1980s, Williamson, Havel and Wiithrich developed 
the idea of extracting the distances between pairs of hydrogen 
atoms in a protein, using nuclear magnetic resonance (NMR). 
The extracted distances were then used to reconstruct 3D 
shapes of molecules^ [5]. The NMR spectrometer (together 
with some post-processing) outputs the distances between 
pairs of atoms in a large molecule. The distances are not 
specified for all atom pairs, and they are uncertain—given 
only up to an interval. This setup lends itself naturally to EDM 
treatment; for example, it can be directly addressed using MDS 
[21]. Indeed, the crystallography community also contributed 
a large number of important results on distance geometry. In a 
different biochemical application, comparing distance matrices 
yields efficient algorithms for comparing proteins from their 
3D structure [22]. 

In machine learning, one can learn manifolds by finding 
an EDM with a low embedding dimension that preserves the 
geometry of local neighborhoods. Weinberger and Saul use 
it to learn image manifolds [9]. Other examples of using 
Euclidean distance geometry in machine learning are results 
by Tenenbaum, De Silva and Langford [23] on image under¬ 
standing and handwriting recognition, Jain and Saul [24] on 
speech and music, and Demaine and et al. [25] on music and 
musical rhythms. 

With the increased interest in sensor networks, several 
EDM-based approaches were proposed for sensor localization 
[2], [3], [4], [20]. Connections between EDMs, multilateration 
and semidefinite programming are expounded in depth in [26], 
especially in the context of sensor network localization. 

Position calibration in ad-hoc microphone arrays is often 
done with sources at unknown locations, such as handclaps, 
fingersnaps or randomly placed loudspeakers [27], [28], [29]. 
This gives us distances (possibly up to an offset time) between 
the microphones and the sources and leads to the problem of 
multi-dimensional unfolding [30]. 

All of the above applications work with labeled distance 
data. In certain TOA- based applications one loses the labels— 
the correct permutation of the distances is no longer known. 
This arises in reconstructing the geometry of a room from 
echoes [6]. Another example of unlabeled distances is in sparse 
phase retrieval, where the distances between the unknown 
non-zero lags in a signal are revealed in its autocorrelation 
function [7]. Recently, motivated by problems in crystallog- 

^Wtithi'ich received the Nobel Prize for chemistry in 2002. 
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raphy, Gujarahati and co-authors published an algorithm for 
reconstruction of Euclidean networks from unlabeled distance 
data [31], 

B. Our Mission 

We were motivated to write this tutorial after realizing that 
EDMs are not common knowledge in the signal processing 
community, perhaps for the lack of a compact introductory 
text. This is effectively illustrated by the anecdote that, not 
long before writing this article, one of the authors had to 
add the (rather fundamental) rank property to the Wikipedia 
page on EDMs."^ In a compact tutorial we do not attempt 
to be exhaustive; much more thorough literature reviews are 
available in longer exposes on EDMs and distance geometry 
[10], [32], [33]. Unlike these works that take the most general 
approach through graph realizations, we opt to show simple 
cases through examples, and to explain and spell out a set of 
basic algorithms that anyone can use immediately. Two big 
topics that we discuss are not commonly treated in the EDM 
literature; localization from unlabeled distances, and multidi¬ 
mensional unfolding (applied to microphone localization). On 
the other hand, we choose to not explicitly discuss the sensor 
network localization (SNL) problem, as the relevant literature 
is abundant. 

Implementations of all the algorithms are available online.^ 
Our hope is that this will provide a good entry point for those 
wishing to learn much more, and inspire new approaches to 
old problems. 

TABLE I 

Summary of notation 


such an expression to provide interesting structural insights. 
We will define notation as it becomes necessary—a summary 
is given in Table I. 

Consider a collection of n points in a d-dimensional Eu¬ 
clidean space, ascribed to the columns of matrix X G 
X = \xi, X2, ■ ■ ■ , Xn]^ Xi G Then the squared distance 
between Xi and Xj is given as 

^ij II ^j\\ J 

where || • || denotes the Euclidean norm. Expanding the norm 
yields 

dij = {xi - XjY{xi - Xj) = xjXi - 2xJXj + xjXj. (2) 
Erom here, we can read out the matrix equation for D = [dij], 
edm(X) = 1 diag(X^X)^ - 2X^X + dmg{X^X)l^, 

(3) 

where 1 denotes the column vector of all ones and diag(A) 
is a column vector of the diagonal entries of A. We see that 
edm(X) is in fact a function of X^X. Eor later reference, it 
is convenient to define an operator JC{G) similar to edm(X), 
that operates directly on the Gram matrix G = X^X, 

IC{G) = diag(G) -2G+ 1 diag(G)^. (4) 

The EDM assembly formula (3) or (4) reveals an important 
property; Because the rank of X is at most d (it has d rows), 
then the rank of X^X is also at most d. The remaining two 
summands in (3) have rank one. By rank inequalities, rank of 
a sum of matrices cannot exceed the sum of the ranks of the 
summands. With this observation, we proved one of the most 
notable facts about EDMs; 


Symbol Meaning 


n 

d 

aij 

D 

edm(X) 
edm(X, Y) 

K{G) 

J 

Aiv 

W 


Number of points (columns) in X = [a;i,..., x„] 
Dimensionality of the Euclidean space 
Element of a matrix A on the ith row and the jth column 
A Euclidean distance matrix 

Euclidean distance matrix created from columns in X 
Matrix containing the squared distances between the 
columns of X and Y 

Euclidean distance matrix created from the Gram matrix G 

Geometric centering matrix 

Restriction of A to non-zero entries in W 

Mask matrix, with ones for observed entries 

Set of real symmetric positive semidefinite matrices in 

I^nXn 


afFdim(X) 

AoB 

^ij 

ei 


Affine dimension of the points listed in X 
Hadamard (entrywise) product of A and B 
Noise corrupting the {i,j) distance 
jth vector of the canonical basis 
Frobenius norm of A, 


II. Erom Points to EDMs and Back 

The principal EDM-related task is to reconstruct the original 
point set. This task is an inverse problem to the simpler 
forward problem of finding the EDM given the points. Thus 
it is desirable to have an analytic expression for the EDM in 
terms of the point matrix. Beyond convenience, we can expect 

^We are working on improving that page substantially. 

^http://lcav.epfl.ch/ivan.dokmanic 


Theorem 1 (Rank of EDMs). Rank of an EDM corresponding 
to points in is at most d + 2. 

This is a powerful theorem; it states that the rank of an 
EDM is independent of the number of points that generate it. 
In many applications, d is three or less, while n can be in the 
thousands. According to Theorem 1, rank of such practical 
matrices is at most five. The proof of this theorem is simple, 
but to appreciate that the property is not obvious, you may try 
to compute the rank of the matrix of non-squared distances. 

What really matters in Theorem 1 is the affine dimension of 
the point set—the dimension of the smallest affine subspace 
that contains the points, denoted by affdim(X). Eor example, 
if the points lie on a plane (but not on a line or a circle) in 
rank of the corresponding EDM is four, not five. This will 
be clear from a different perspective in the next subsection, as 
any affine subspace is just a translation of a linear subspace. 
An illustration for a ID subspace of is provided in Eig. 
2; Subtracting any point in the affine subspace from all its 
points translates it to the parallel linear subspace that contains 
the zero vector. 

A. Essential Uniqueness 

When solving an inverse problem, we need to understand 
what is recoverable and what is forever lost in the forward 
problem. Representing sets of points by distances usually 
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Fig. 2. Illustration of the relationship between an affine subspace and its 
parallel linear subspace. The points X = [tci,..., 334 ] live in an affine 
subspace—a line in that does not contain the origin. In (A), the vector 3:1 
is subtracted from all the points, and the new point list is X' = [ 0,352 — 
3:1 , 353 — 351,354 — xi]. While the columns of X span R^, the columns of 
X' only span the ID subspace of —the line through the origin. In (B), we 
subtract a different vector from all points: the centroid ^X1. The translated 
vectors X" = [3:",..., 3:"] again span the same ID subspace. 



Fig. 3. Illustration of a rigid transformation in 2D. Here the points set is 
transformed as RX + bl^. Rotation matrix R= [_5 q], corresponds to a 
counterclockwise rotation of 90°. The translation vector is 6 = [3, 1]^. The 
shape is drawn for visual reference. 

increases the size of the representation. For most interesting n 
and d, the number of pairwise distances is larger than the size 
of the coordinate description, ( 2 ) > nd, so an EDM holds 
more scalars than the list of point coordinates. Nevertheless, 
some information is lost in this encoding, namely the infor¬ 
mation about the absolute position and orientation of the point 
set. Intuitively, it is clear that rigid transformations (including 
reflections) do not change distances between the fixed points 
in a point set. This intuitive fact is easily deduced from the 
EDM assembly formula (3). We have seen in (3) and (4) that 
edm(X) is in fact a function of the Gram matrix X. 

This makes it easy to show algebraically that rotations and 
reflections do not alter the distances. Any rotation/reflection 
can be represented by an orthogonal matrix Q G acting 

on the points Xi. Thus for the rotated point set Xr — QX 
we can write 

XjXr = {QX)^{QX) = X^Q^QX = X^X, (5) 

where we invoked the orthogonality of the rotation/reflection 
matrix, Q = I. 

Translation by a vector 6 S can be expressed as 

Xt=X + bl^. (6) 


Using diag(X7Xt) = diag(X^X) + 2X^b + \\bf 1, one 
can directly verify that this transformation leaves (3) intact. In 
summary, 

edm(QX) = edm(X + bl^) = edm(X). (7) 

The consequence of this invariance is that we will never 
be able to reconstruct the absolute orientation of the point 
set using only the distances, and the corresponding degrees 
of freedom will be chosen freely. Different reconstruction 
procedures will lead to different realizations of the point set, 
all of them being rigid transformations of each other. Fig. 3 
illustrates a point set under a rigid transformation. It is clear 
that the distances between the points are the same for all three 
shapes. 

B. Reconstructing the Point Set From Distances 

The EDM equation (3) hints at a procedure to compute 
the point set starting from the distance matrix. Consider the 
following choice: let the first point Xi be at the origin. Then 
the first column of D contains the squared norms of the point 
vectors, 

dii = \\xi - XiW^ = \\xi - of' = . ( 8 ) 

Consequently, we can immediately construct the term 
1 diag(X^X) and its transpose in (3), as the diagonal of 
X^X contains exactly the norms squared ||a;i||^. Concretely, 

ld\ag{X^X) = Idl, (9) 

where di = Dei is the first column of D. We thus obtain 
the Gram matrix from (3) as 

G = X^X =-]^{D - IdJ-dil^). ( 10 ) 

The point set can be found by an EVD, G = UAU^, 
where A = diag(Ai,..., A„) with all eigenvalues Xi non¬ 
negative, and U orthonormal, as G is a symmetric posi¬ 
tive semideflnite matrix. Throughout the paper we assume 
that the eigenvalues are sorted in the order of decreasing 
magnitude, |Ai| > IA 2 I > ••• > |A„|. We can now set 
X = [diag(5/Ar, fXd), Odx(n-d)] Note that we 
could have simply taken A^^^U^ as the reconstructed point 
set, but if the Gram matrix really describes a d-dimensional 
point set, the trailing eigenvalues will be zeroes, so we choose 
to truncate the corresponding rows. 

It is straightforward to verify that the reconstructed point 
set X generates the original EDM, D = edm(X); as we 
have learned, X and X are related by a rigid transformation. 
The described procedure is called the classical MDS, with a 
particular choice of the coordinate system: a:i is fixed at the 
origin. 

In (10) we subtract a structured rank-2 matrix {IdJ -f 
dil^) from D. A more systematic approach to the classical 
MDS is to use a generalization of ( 10 ) by Gower [ 16 ]. Any 
such subtraction that makes the right hand side of ( 10 ) positive 
semideflnite (PSD), i.e., that makes G a Gram matrix, can also 
be modeled by multiplying D from both sides by a particular 
matrix. This is substantiated in the following result. 










5 


Algorithm 1 Classical MDS 
1: function ClassicalMDS(D, (f) 

2: J ^ I — ^11^ \> Geometric centering matrix 

3: G ^ k JDJ > Compute the Gram matrix 

4: 17,[A,]r=i^EVD(G) 

5: return [diag(yA];, ^Ad), Odx(n-d)\U^ 

6: end function 


Theorem 2 (Gower [16]). D is an EDM if and only if 

ls^)D{I-sl^) (11) 

is PSD for any s such that 1 = 1 and D f 0. 

In fact, if (11) is PSD for one such s, then it is PSD for all 
of them. In particular, define the geometric centering matrix 
as 

J = (12) 

n 

Then ~^JDJ being positive semidefinite is equivalent to D 
being an EDM. Different choices of s correspond to different 
translations of the point set. 

The classical MDS algorithm with the geometric centering 
matrix is spelled out in Algorithm 1. Whereas so far we 
have assumed that the distance measurements are noiseless. 
Algorithm 1 can handle noisy distances too, as it discards all 
but the d largest eigenvalues. 

It is straightforward to verify that (10) corresponds to s = 
ei. Think about what this means in terms of the point set: 
Xei selects the first point in the list, Xi. Then Xq = X{I — 
eil^) translates the points so that Xi is translated to the 
origin. Multiplying the definition (3) from the right by (J — 
eil^) and from the left by {I — lej) will annihilate the two 
rank-1 matrices, diag(G) 1 ^ and 1 diag(G)^. We see that the 
remaining term has the form —2Xq Xq, and the reconstructed 
point set will have the first point at the origin! 

On the other hand, setting s = -1 places the centroid of 
the point set at the origin of the coordinate system. For this 
reason, the matrix J = I — ^11^ is called the centering 
matrix. To better understand why, consider how we normally 
center a set of points given in X. 

First, we compute the centroid as the mean of all the points 

1 n ^ 

Xc= -y^ Xi= -XI. (13) 

i=l 

Second, we subtract this vector from all the points in the set, 

X^ = X-x^l^ =X--X11^ = X{I--11^). (14) 
n n 

In complete analogy with the reasoning for s = ei, we can see 
that the reconstructed point set will be centered at the origin. 

C. Orthogonal Procrustes Problem 

Since the absolute position and orientation of the points 
are lost when going over to distances, we need a method to 
align the reconstructed point set with a set of anchors —^points 
whose coordinates are fixed and known. 


This can be achieved in two steps, sometimes called Pro¬ 
crustes analysis. Ascribe the anchors to the columns of Y, 
and suppose that we want to align the point set X with the 
columns of Y . Let Xa denote the submatrix (a selection of 
columns) of X that should be aligned with the anchors. We 
note that the number of anchors—columns in Xa —is typically 
small compared with the total number of points—columns in 
X. 

In the first step, we remove the means j/c and Xa^c from 
matrices Y and Xa, obtaining the matrices Y and Xa. In 
the second step, termed orthogonal Procrustes analysis, we 
are searching for the rotation and reflection that best maps 
Xa onto Y, 

R= argmin ||QXa —(15) 

Q:QQ^=I 

The Frobenius norm || • ||p, is simply the f^-norm of the matrix 
entries, ||A||^ ~ trace(A^A). 

The solution to (15)— found by Schonemann in his PhD the¬ 
sis [34]— is given by the singular value decomposition (SVD). 
Let XaY — then we can continue computing (15) 

as follows 

R= argmin ||(5Xa||^-f || — trace(Y^(5Xa) 

Q.QQ^^I 

= argmax trace(QS), (16) 

Q.QQ^^I 

where Q V^QU, and we used the orthogonal invariance 
of the Frobenius norm and the cyclic invariance of the trace. 
The last trace expression in (16) is equal to X]r=r ^idu- Noting 
that Q is also an orthogonal matrix, its diagonal entries cannot 
exceed 1. Therefore, the maximum is achieved when qa = 1 
for all i, meaning that the optimal Q is an identity matrix. It 
follows that R = VU^. 

Once the optimal rigid transformation has been found, the 
alignment can be applied to the entire point set as 

R{X - Xa^cl^) . (17) 


D. Counting the Degrees of Freedom 

It is interesting to count how many degrees of freedom there 
are in different EDM related objects. Clearly, for n points in 
we have 

#x = n X d (18) 

degrees of freedom: If we describe the point set by the list 
of coordinates, the size of the description matches the number 
of degrees of freedom. Going from the points to the EDM 
(usually) increases the description size to ^n{n — 1), as the 
EDM lists the distances between all the pairs of points. By 
Theorem 1 we know that the EDM has rank at most d + 2. 

Let us imagine for a moment that we do not know any 
other EDM-specific properties of our matrix, except that it 
is symmetric, positive, zero-diagonal (or hollow), and that 
it has rank d + 2. The purpose of this exercise is to count 
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the degrees of freedom associated with such a matrix, and 
to see if their number matches the intrinsic number of the 
degrees of freedom of the point set, If it did> then 

these properties would completely characterize an EDM. We 
can already anticipate from Theorem 2 that we need more 
properties: a certain matrix related to the EDM—as given in 
(11) —must be PSD. Still, we want to see how many degrees 
of freedom we miss. 

We can do the counting by looking at the EVD of a 
symmetric matrix, D = UAU^ . The diagonal matrix A is 
specihed by d + 2 degrees of freedom, because D has rank 
d +2. The hrst eigenvector of length n takes up n — 1 degrees 
of freedom due to the normalization; the second one takes up 
n — 2, as it is in addition orthogonal to the hrst one; for the 
last eigenvector, number (d + 2), we need n — (d + 2) degrees 
of freedom. We do not need to count the other eigenvectors, 
because they correspond to zero eigenvalues. The total number 
is then 


^DOF = (d + 2) + (n — 1) + • • • + [n — (d + 2)] — n 

Eigenvalues Eigenvectors Hollowness 

. . (d + 1) X (d + 2) 

= n X (d + 1) - i^i 


For large n and fixed d, it follows that 


^DOF d 

#x d 


(19) 


Therefore, even though the rank property is useful and we will 
show efficient algorithms that exploit it, it is still not a tight 
property (symmetry and hollowness included). Eor d = 3, 
the ratio (19) is |, so loosely speaking, the rank property 
has 30% determining scalars too many, which we need to 
set consistently. Put differently, we need 30% more data in 
order to exploit the rank property than we need to exploit the 
full EDM structure. Again loosely phrased, we can assert that 
for the same amount of data, the algorithms perform at least 
r:! 30% worse if we only exploit the rank property, without 
EDMness. 

The one-third gap accounts for various geometrical con¬ 
straints that must be satished. The redundancy in the EDM 
representation is what makes denoising and completion algo¬ 
rithms possible, and thinking in terms of degrees of freedom 
gives us a fundamental understanding of what is achievable. 
Interestingly, the above discussion suggests that for large n 
and large d = o(n), little is lost by only considering rank. 

Einally, in the above discussion, for the sake of simplicity 
we ignored the degrees of freedom related to absolute orien¬ 
tation. These degrees of freedom, not present in the EDM, do 
not affect the large-n behavior. 


E. Summary 

Let us summarize what we have achieved in this section: 

• We explained how to algebraically construct an EDM 
given the list of point coordinates. 


• We discussed the essential uniqueness of the point set; 
information about the absolute orientation of the points 
is irretrievably lost when transitioning from points to an 
EDM, 

• We explained classical MDS—a simple eigenvalue- 
decomposition-based algorithm (Algorithm 1) for re¬ 
constructing the original points—along with discussing 
parameter choices that lead to different centroids in 
reconstruction, 

• Degrees-of-freedom provide insight into scaling behavior. 
We showed that the rank property is pretty good, but there 
is more to it than just rank. 


III. EDMs AS A Practical Tool 


We rarely have a perfect EDM. Not only are the entries 
of the measured matrix plagued by errors, but often we can 
measure just a subset. There are various sources of error 
in distance measurements: we already know that in NMR 
spectroscopy, instead of exact distances we get intervals. 
Measuring distance using received powers or TOAs is subject 
to noise, sampling errors and model mismatch. 

Missing entries arise because of the limited radio range, 
or because of the nature of the spectrometer. Sometimes the 
nodes in the problem at hand are asymmetric by dehnition; in 
microphone calibration we have two types: microphones and 
calibration sources. This results in a particular block structure 
of the missing entries (we will come back to this later, but 
you can fast-forward to Eig. 5 for an illustration). 

It is convenient to have a single statement for both EDM ap¬ 
proximation and EDM completion, as the algorithms described 
in this section handle them at once. 

Problem 1. Let D = edm(X). We are given a noisy 
observation of the distances between p < ( 2 ) pairs of points 
from X. That is, we have a noisy measurement of 2p entries 
in D, 

dij = dij Cij , ( 20 ) 

for {i,j) € E, where E is some index set, and absorbs 
all errors. The goal is to reconstruct the point set X in the 
given embedding dimension, so that the entries of ed'm(X) 
are close in some metric to the observed entries dij. 

To concisely write down completion problems, we dehne 
the mask matrix IV as follows. 


def 

W,j = 


1, (i,j) € E 
0, otherwise. 


( 21 ) 


This matrix then selects elements of an EDM through a 
Hadamard (entrywise) product. Eor example, to compute the 
norm of the difference between the observed entries in A and 
B, we write ||W o(A —B)||. Eurthermore, we dehne the 
indexing Avk to mean the restriction of A to those entries 
where W is non-zero. The meaning of Bw ^ Aw is that 
we assign the observed part of A to the observed part of B. 
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Algorithm 2 Alternating Rank-Based EDM Completion 

1: function RankCompleteEDM(VE, £), d) 

2: Dw -i— Dw \> Initialize observed entries 

3: DiiT -w ^ M Initialize unobserved entries 

4: repeat 

5: EVThreshold(iI>, d -I- 2) 

6: Dw -i— £)vv i> Enforce known entries 

7: Di -i— 0 > Set the diagonal to zero 

8: Z) -i— (-D)+ > Zero the negative entries 

9: until Convergence or Maxiter 

10 : return D 

11 : end function 

12: function EVTHRESHOLD(D,r) 

13: U, ^ EVD(Z>) 

14: S ^ diag (Ai, ■. ■, Xr, 0, ■ - > Oj 

15: D 3- ri—r times 

16: return D 

17: end function 


A. Exploiting the Rank Property 

Perhaps the most notable fact about EDMs is the rank prop¬ 
erty established in Theorem 1 : The rank of an EDM for points 
living in is at most d-|-2. This leads to conceptually simple 
algorithms for EDM completion and denoising. Interestingly, 
these algorithms exploit only the rank of the EDM. There is 
no explicit Euclidean geometry involved, at least not before 
reconstructing the point set. 

We have two pieces of information; a subset of potentially 
noisy distances, and the desired embedding dimension of the 
point configuration. The latter implies the rank property of the 
EDM that we aim to exploit. We may try to alternate between 
enforcing these two properties, and hope that the algorithm 
produces a sequence of matrices that converges to an EDM. If 
it does, we have a solution. Alternatively, it may happen that 
we converge to a matrix with the correct rank that is not an 
EDM, or that the algorithm never converges. The pseudocode 
is listed in Algorithm 2. 

A different, more powerful approach is to leverage al¬ 
gorithms for low rank matrix completion developed by the 
compressed sensing community. Eor example, OptSpace [35] 
is an algorithm for recovering a low-rank matrix from noisy, 
incomplete data. Let us take a look at how OptSpace works. 
Denote by M G the rank-r matrix that we seek 

to recover, by Z G measurement noise, and by 

W G the mask corresponding to the measured entries; 

for simplicity we choose m < n. The measured noisy and 
incomplete matrix is then given as 

M = W o{M + Z). (22) 


Effectively, this sets the missing (non-observed) entries of the 
matrix to zero. OptSpace aims to minimize the following cost 
function. 


F{A,S,B) = l Wo{M-ASB^)^, (23) 

Z F 


Algorithm 3 OptSpace [35] 

1 : function OPTSP^E(M,r) 

2: M G- Trim(7Vi’) 

3: ^ SVD(a-iM) ^ 

4: Aq ^ Eirst r columns of A 

5: Bq g- Eirst r columns of B 

6: So ■<— cLigmmF(Ao, S, Bq) [> Eq. (23) 

S^Rr-Xr 

7: A^B ^ argmin_F(A, S'o, B) > See the note below 

B—I 

8 : return ASqB^ 

9: end function 

> Line 7: gradient descent starting at Aq, Bq 


where S' S A G and B G such that 

A^A — B^B — I. Note that S need not be diagonal. 

The cost function (23) is not convex, and minimizing it is a 
priori difficult [36] due to many local minima. Nevertheless, 
Keshavan, Montanari and Oh [35] show that using the gradient 
descent method to solve (23) yields the global optimum 
with high probability, provided that the descent is correctly 
initialize _ 

Let M = be the SVD of M. Then 

'—' '—' def 

we define the scaled rank-r projection of M as Mr = 
a~^ 'Fg=i ■ The fraction of observed entries is denoted 

by a, so that the scaling factor conyjensates the smaller 
average magnitudejrf the entries in M in comparison with 
M. The SVD of Mr is then used to initialize the gradient 
descent, as detailed in Algorithm 3. 

Two additional remarks are due in the description of 
OptSpace. Eirst, it can be shown that the performance is 
improved by zeroing the over-represented rows and columns. 
A row (resp. column) is over-represented if it contains more 
than twice the average number of observed entries per row 
{resp. column). These heavy rows and columns bias the cor¬ 
responding singular vectors and singular values, so (perhaps 
surprisingly) it is better to throw them away. We call this step 
“Trim” in Algorithm 3. 

Second, the minimization of (23) does not have to be 
performed for all variables at once. In [35], the authors first 
solve the easier, convex minimization for S, and then with the 
optimizer S fixed, they find the matrices A and B using the 
gradient descent. These steps correspond to lines 6 and 7 of 
Algorithm 3. Eor an application of OptSpace in calibration 
of ultrasound measurement rigs, see the “Calibration in 
Ultrasound Tomography” box. 

B. Multidimensional Scaling 

Multidimensional scaling refers to a group of techniques 
that, given a set of noisy distances, find the best fitting point 
conformation. It was originally proposed in psychometrics 
[15], [8] to visualize the (dis-)similarities between objects. 
Initially, MDS was defined as the problem of representing 
distance data, but now the term is commonly used to refer 
to methods for solving the problem [39]. 

Various cost functions were proposed for solving MDS. In 
Section II-B, we already encountered one method: the classical 
















Calibration in Ultrasound Tomography 


The rank property of EDMs, introduced in Theorem 1 can be 
leveraged in calibration of ultrasound tomography devices. An 
example device for diagnosing breast cancer is a circular ring w/ith 
thousands of ultrasound transducers, placed around the breast 
[37], The setup is shown in Fig. 4A. 

Due to manufacturing errors, the sensors are not located on 
a perfect circle. This uncertainty in the positions of the sen¬ 
sors negatively affects the algorithms for imaging the breast. 
Fortunately, we can use the measured distances between the 
sensors to calibrate their relative positions. We can estimate the 
distances by measuring the times-of-flight (TOFs) between pairs 
of transducers in a homogeneous environment, e.g. in water. 

We cannot estimate the distances between all pairs of sensors 
because the sensors have limited beam widths (it is hard to 
manufacture omni-directional ultrasonic sensors). Therefore, the 
distances between neighboring sensors are unknown, contrary to 
typical SNL scenarios where only the distances between nearby 
nodes can be measured. Moreover, the distances are noisy and 
some of them are unreliably estimated. This yields a noisy and 
incomplete EDM, whose structure is illustrated in Figure 46. 
Assuming that the sensors lie in the same plane, the original 
EDM produced by them would have a rank less than five. We can 
use the rank property and a low-rank matrix completion method, 
such as OptSpace (Algorithm 3), to complete and denoise the 
measured matrix [38]. Then we can use the classical MDS in 
Algorithm 1 to estimate the relative locations of the ultrasound 


sensors. 

For reasons mentioned above, SNL-specific algorithms are sub- 
optimal when applied to ultrasound calibration. An algorithm 
based on the rank property effectively solves the problem, and 
enables one to derive upper bounds on the performance error 
calibration mechanism, with respect to the number of sensors 
and the measurement noise. The authors in [38] show that the 
error vanishes as the number of sensors increases. 



Fig. 4. (A) Ultrasound transducers lie on an approximately circular 

ring. The ring surrounds the breast and after each transducer fires 
an ultrasonic signal, the sound speed distribution of the breast is 
estimated. A precise knowledge of the sensor locations is needed to 
have an accurate reconstruction of the enclosed medium. (B) Because 
of the limited beam width of the transducers, noise and imperfect TOF 
estimation methods, the measured EDM is incomplete and noisy. Gray 
areas show missing entries of the matrix. 


MDS. This method minimizes the Frobenius norm of the 
difference between the input matrix and the Gram matrix of 
the points in the target embedding dimension. 

The Gram matrix contains inner products; rather than with 
inner products, it is better to directly work with the distances. 
A typical cost function represents the dissimilarity of the 
observed distances and the distances between the estimated 
point locations. An essential observation is that the feasible set 
for these optimizations is not convex (EDMs with embedding 
dimensions smaller than n — 1 lie on the boundary of a cone 
[20], which is a non-convex set). 

A popular dissimilarity measure is raw stress [40], defined 
as the value of 


minimize 


(iJ)eE ^ 



(24) 


where E defines the set of revealed elements of the distance 
matrix D. The objective function can be concisely written as 
||VFo(^edm(X)-\/B)|| a drawback of this cost function 
is that it is not globally differentiable. Approaches described 
in the literature comprise iterative majorization [41], various 
methods using convex analysis [42] and steepest descent 
methods [43]. 

Another well-known cost function is s-stress. 


minimize 

XgR'^X’* 


(edm(X)y - d.j'j . 


(25) 


Again, we write the objective concisely as ||TF o (edm(X) — 
iD)||^. It was first studied by Takane, Young and De Leeuw 


[44]. Conveniently, the s-stress objective is everywhere differ¬ 
entiable, but at a disadvantage that it favors large over small 
distances. Gaffke and Mathar [45] propose an algorithm to find 
the global minimum of the s-stress function for embedding 
dimension d = n — 1. EDMs with this embedding dimension 
exceptionally constitute a convex set [20], but we are typically 
interested in embedding dimensions much smaller than n. The 
s-stress minimization in (25) is not convex for d < n—1. It was 
analytically shown to have saddle points [46], but interestingly, 
no analytical non-global minimizer has been found [46]. 

Browne proposed a method for computing s-stress based 
on Newton-Raphson root finding [47]. Glunt reports that the 
method by Browne converges to the global minimum of (25) 
in 90% of the test cases in his dataset^ [48]. 

The cost function in (25) is separable across points i 
and across coordinates k, which is convenient for distributed 
implementations. Parhizkar [46] proposed an alternating co¬ 
ordinate descent method that leverages this separability, by 
updating a single coordinate of a particular point at a time. 
The s-stress function restricted to the A:th coordinate of the *th 
point is a fourth-order polynomial, 

f{x; ^ (26) 

where lists the polynomial coefficients for ith point 

and kth coordinate. Eor example, '^ij^ that is, 

four times the number of points connected to point i. Expres¬ 
sions for the remaining coefficients are given in [46]; in the 

®While the experimental setup of Glunt [48] is not detailed, it was 
mentioned that the EDMs were produced randomly. 












9 


Algorithm 4 Alternating Descent [46] 

1: function AlternatingDescent(£), W, d) 

2: X e ^ Xo = 0 > Initialize the point set 

3: repeat 

4: for i S {1, • • • , n} do > Points 

5: for fc G {1, • • • , d} do t> Coordinates 

6: ^ GetQuadricCoeffs(VF, 13, d) 

7 : Xi^k ^ argmin^j, f{x; > Eq. (26) 

8 : end for 

9: end for 

10: until Convergence or Maxiter 

11 : return X 

12 : end function 


pseudocode (Algorithm 4), we assume that these coefficients 
are returned by the function “GetQuadricCoeffs”, given the 
noisy incomplete matrix D, the observation mask W and the 
dimensionality d. The global minimizer of (26) can be found 
analytically by calculating the roots of its derivative (a cubic). 
The process is then repeated over all coordinates k, and points 
i, until convergence. The resulting algorithm is remarkably 
simple, yet empirically converges fast. It naturally lends itself 
to a distributed implementation. We spell it out in Algorithm 
4. 

When applied to a large dataset of random, noiseless and 
complete distance matrices. Algorithm 4 converges to the 
global minimum of (25) in more than 99% of the cases [46]. 


C. Semidefinite Programming 

Recall the characterization of EDMs (11) in Theorem 2. It 
states that D is an EDM if and only if the corresponding 
geometrically centered Gram matrix —^JDJ is positive- 
semidehnite. Thus, it establishes a one-to-one correspondence 
between the cone of EDMs, denoted by EDM”, and the 
intersection of the symmetric positive-semidefinite cone S” 
with the geometrically centered cone S”. The latter is dehned 
as the set of all symmetric matrices whose column sum 
vanishes. 


S” = {G e M”""” I G = G3 = 0} . (27) 

We can use this correspondence to cast EDM completion 
and approximation as semidefinite programs. While (11) de¬ 
scribes an EDM of an n-point configuration in any dimension, 
we are often interested in situations where d n.lt is easy to 
adjust for this case by requiring that the rank of the centered 
Gram matrix be bounded. One can verify that 

£> = edm(X) 1 ^ 

afFdim(X) < dj |rank(J£)J) < d, 

when n > d. That is, EDMs with a particular embedding 
dimension d are completely characterized by the rank and 
dehniteness of JDJ. 


Now we can write the following rank-constrained semidef¬ 
inite program for solving Problem 1, 


minimize 

G 

subject to 


PEo (i3-/C(G)) 
rank(G) < d 


2 


F 


Ge S” nS”. 


(29) 


The second constraint is just a shorthand for writing G ^ 
0, G1 — 0. We note that this is equivalent to MDS with the 
s-stress cost function, thanks to the rank characterization (28). 

Unfortunately, the rank property makes the feasible set in 
(29) non-convex, and solving it exactly becomes difficult. 
This makes sense, as we know that s-stress is not convex. 
Nevertheless, we may relax the hard problem, by simply 
omitting the rank constraint, and hope to obtain a solution 
with the correct dimensionality. 


minimize 

G 

subject to 


Wo[p-lC{G)^ 
Ge S!;:nS”. 


2 

F 


(30) 


We call (30) a semidefinite relaxation (SDR) of the rank- 
constrained program (29). 

The constraint G G S”, or equivalently, G1 — 0, means 
that there are no strictly positive dehnite solutions (G has a 
nullspace, so at least one eigenvalue must be zero). In other 
words, there exist no strictly feasible points [32]. This may 
pose a numerical problem, especially for various interior point 
methods. The idea is then to reduce the size of the Gram matrix 
through an invertible transformation, somehow removing the 
part of it responsible for the nullspace. In what follows, we 
describe how to construct this smaller Gram matrix. 

A different, equivalent way to phrase the multiplicative 
characterization (11) is the following statement: a symmetric 
hollow matrix D is an EDM if and only if it is negative 
semidefinite on {J}^ (on all vectors t such that F1 — 0). 
Let us construct an orthonormal basis for this orthogonal 
complement—a subspace of dimension (n — 1)—and arrange 
it in the columns of matrix V G We demand 


V^l = 0 
V^V = 1. 


(31) 


There are many possible choices for V, but all of them obey 
that =1—-lF = J. The following choice is given 

in [2], 



p 

p 

P 


1 + g q 

q 

V = 

q l + q ■■■ 

q 


q 

q 

l + q_ 


where p = — l/(n -f ^Jn) and q = —lly/n. 

With the help of the matrix V, we can now construct the 
sought Gramian with reduced dimensions. Eor an EDM D G 

j^nxn 

g{D) = -]^V^DV 


(33) 
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Algorithm 5 Semidefinite Relaxation (Matlab/CVX) 

1 

function [EDM, X] = 

sdr_complete. 

_edm(D, W, lambda) 

3 

n = 

size(D, 1 ); 



4 

X = 

-1/ {n + sqrt(n) ) 

; 


5 

y = 

-1/sqrt(n) ; 



6 

V = 

[y*ones ( 1, n-1 ); 

x*ones (n-1) 

+ eye (n-1 )]; 

7 

e = 

ones (n, 1 ); 



9 

cvx_ 

_begin sdp 



10 


variable G(n-1, 

n-1) symmetric; 

11 


B = v*G*V' ; 



12 


E = diag(B)*e ' + 

e*diag (B) ' 

- 2*B; 

13 


maximize trace(G) ... 


14 


- lambda 

* norm ( W .* 

(E - D) , ' fro ' ) ; 

15 


subject to 



16 


G >= 0; 



17 

cvx_ 

_end 



19 

[U, 

S, V] = svd(B) ; 



20 

EDM 

= diag(B)*e ' + e 

*diag(B) ' - 

2*B; 

21 

X = 

sqrt (S)*V'; 




is an (n — 1) x (n — 1) PSD matrix. This can be verified by 
substituting (33) in (4). Additionally, we have that 

K.{Vg{D)V^) = D. (34) 

Indeed, H i— JCiVHV^) is an invertible mapping from 
to EDM” whose inverse is exactly Q. Using these nota¬ 
tions we can write down an equivalent optimization program 
that is numerically more stable than (30) [2]: 


Microphones Acoustic events 



D = 






Fig. 5. Microphone calibration as an example of MDU. We can measure 
only the propagation times from acoustic sources at unknown locations, to 
microphones at unknown locations. The corresponding revealed part of the 
EDM has a pai'ticular off-diagonal structure, leading to a special case of EDM 
completion. 


Noting that trace(ff) = trace(G) because trace(J£)J) = 
trace(V^£)V), we write the following SDR, 


maximize 

H 

subject to 


trace(iT) — A 

H e 


Wo{D-K,{VHV^)) 

F 

(36) 


Here we opted to include the data fidelity term in the La- 
grangian form, as proposed by Biswas [49], but it could also be 
moved to constraints. Finally, in all of the above relaxations, 
it is straightforward to include upper and lower bounds on 
the distances. Because the bounds are linear constraints, the 
resulting programs remain convex; this is particularly useful 
in the molecular conformation problem. A Matlab/CVX [50], 
[51] implementation of the SDR (36) is given in Algorithm 5. 


minimize 

H 

subject to 


W o[d-1C{VHV^)) 

H G 


2 


F 


(35) D. Multidimensional Unfolding: A Special Case of Comple¬ 
tion 


On the one hand, with the above transformation the constraint 
G1 = 0 became implicit in the objective, as VHV^ 1 = 0 
by (31); on the other hand, the feasible set is now the full 
semidefinite cone S"“^. 

Still, as Krislock & Wolkowicz mention [32], by omitting 
the rank constraint we allow the points to move about in a 
larger space, so we may end up with a higher-dimensional 
solution even if there is a completion in dimension d. 

There exist various heuristics for promoting lower rank. One 
such heuristic involves the trace norm—the convex envelope 
of rank. The trace or nuclear norm is studied extensively by 
the compressed sensing community. In contrast to the common 
wisdom in compressed sensing, the trick here is to maximize 
the trace norm, not to minimize it. The mechanics are as 
follows: maximizing the sum of squared distances between 
the points will stretch the configuration as much as possible, 
subject to available constraints. But stretching favors smaller 
affine dimensions (imagine pulling out a roll of paper, or 
stretching a bent string). Maximizing the sum of squared 
distances can be rewritten as maximizing the sum of norms 
in a centered point configuration—^but that is exactly the 
trace of the Gram matrix G = —^JDJ [9]. This idea has 
been successfully put to work by Weinberger and Saul [9] in 
manifold learning, and by Biswas et al. in SNL [49]. 


Imagine that we partition the point set into two subsets, 
and that we can measure the distances between the points 
belonging to different subsets, but not between the points in 
the same subset. Metric multidimensional unfolding (MDU) 
[30] refers to this special case of EDM completion. 

MDU is relevant for position calibration of ad-hoc sensor 
networks, in particular of microphones. Consider an ad-hoc ar¬ 
ray of m microphones at unknown locations. We can measure 
the distances to k point sources, also at unknown locations, for 
example by emitting a pulse (we assume that the sources and 
the microphones are synchronized). We can always permute 
the points so that the matrix assumes the structure shown in 
Fig. 5, with the unknown entries in two diagonal blocks. This 
is a standard scenario described for example in [27]. 

One of the early approaches to metric MDU is that of 
Schonemann [30]. We go through the steps of the algorithm, 
and then explain how to solve the problem using the EDM 
toolbox. The goal is to make a comparison, and emphasize 
the universality and simplicity of the introduced tools. 

Denote by i? = [ri, ..., r^j the unknown microphone 
locations, and by S' = [si, ..., Sfc] the unknown source 
locations. The distance between the jth microphone and jth 
source is 


~ IIG Sjll , 


(37) 
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Fig. 6 . Comparison of different algorithms applied to completing an EDM 
with random deletions. For every number of deletions, we generated 2000 
realizations of 20 points uniformly at random in a unit square. Distances to 
delete were chosen uniformly at random among the resulting (^ 2 ^) = 190 
pairs; 20 deletions correspond to 5 ^ 10% of the number of distance pairs and 
to 5% of the number of matrix entries; 150 deletions correspond to 5 ^ 80% 
of the distance pairs and to ^ 38% of the number of matrix entries. Success 
was declared if the Frobenius norm of the error between the estimated matrix 
and the true EDM was less than 1% of the Frobenius norm of the true EDM. 



Number of acoustic events 


Fig. 7. Comparison of different algorithms applied to multidimensional 
unfolding with varying number of acoustic events k. For every number of 
acoustic events, we generated 3000 realizations of m = 20 microphone 
locations uniformly at random in a unit cube. Percentage of missing matrix 
entries is given as + 'm?)/{k + m)^, so that the ticks on the abscissa 
correspond to [68,56,51,50,51,52]% (non-monotonic in k with the mini¬ 
mum for k = m = 20). Success was declared if the Frobenius norm of the 
en'or between the estimated matrix and the true EDM was less than 1% of 
the Frobenius norm of the true EDM. 


SO that in analogy with (3) we have 

A = edm{R,S) = diag{R^ R)l^-2R^ S+1 diag{S^ S), 

(38) 

where we overloaded the edm operator in a natural way. We 
use A to avoid confusion with the standard Euclidean D. 
Consider now two geometric centering matrices of sizes m 
and k, denoted Jm and Similarly to (14), we have 

RJm = R-r,l^. SJu = S-s^l^. (39) 

This means that 

JmAJk = R^S = G (40) 

is a matrix of inner products between vectors and Sj. 
We used tildes to differentiate this from real inner products 
betwen and Sj, because in (40), the points in R and S are 
referenced to different coordinate systems. The centroids 
and Sc generally do not coincide. There are different ways to 
decompose G into a product of two full rank matrices, call 
them A and B, 

G = A^B. (41) 

We could for example use the SVD, G = , and set 

A^ = U and B = . Any two such decompositions are 

linked by some invertible transformation T G 

G = A^B = R^T-^TS. (42) 

We can now write down the conversion rule from what we 
can measure to what we can compute, 

R = T^A + rcl^ 

, X X (43) 

S = {T-y B + Scl^, 

where A and B can be computed according to (41). Because 
we cannot reconstruct the absolute position of the point set, 
we can arbitrarily set = 0, and s^ = aei. Recapitulating, 
we have that 

A = edur{T^A, {T-^YB + aeil^) , (44) 


and the problem is reduced to computing T and a so that 
(44) hold, or in other words, so that the right hand side be 
consistent with the data A. We reduced MDU to a relatively 
small problem: in 3D, we need to compute only ten scalars. 
Schonemann [30] gives an algebraic method to find these 
parameters, and mentions the possibility of least squares, while 
Crocco, Bue and Murino [27] propose a different approach 
using non-linear least squares. 

This procedure seems quite convoluted. Rather, we see 
MDU as a special case of matrix completion, with the structure 
illustrated in Fig. 5. 

More concretely, represent the microphones and the sources 
by a set of n = k-\-m points, ascribed to the columns of matrix 
JV = [i? S']. Then edm(X) has a special structure as seen in 
Fig. 5, 


edm(X) 


edm(i?) edm(il, S) 
edm(S, ii) edm(S) 


We define the mask matrix for MDU as 


(45) 


XT/ def 

*Vmdu — 


0 rnxm 
^ kxni 


^ mxk 
Okxk 


(46) 


With this matrix, we can simply invoke the SDR in Algorithm 
5. We could also use Algorithm 2, or Algorithm 4. Perfor¬ 
mance of different algorithms is compared in Section III-E. 

It is worth mentioning that SNL specific algorithms that 
exploit the particular graph induced by limited range commu¬ 
nication do not perform well on MDU. This is because the 
structure of the missing entries in MDU is in a certain sense 
opposite to the one of SNL. 


E. Performance Comparison of Algorithms 

We compare the described algorithms in two different EDM 
completion settings. In the first experiment (Figs. 6 and 8), the 
entries to delete are chosen uniformly at random. The second 
experiment (Figs. 7 and 9) tests performance in MDU, where 
the non-observed entries are highly structured. In Figs. 6 and 
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Number of deletions 


Number of acoustic events 


Fig. 8. Comparison of different algorithms applied to completing an EDM 
with random deletions and noisy distances. For every number of deletions, we 
generated 1000 realizations of 20 points uniformly at random in a unit square. 
In addition to the number of deletions, we varied the amount of jitter added 
to the distances. Jitter was drawn from a centered uniform distribution, with 
the level increasing in the direction of the arrow, from W[0, 0] (no jitter) for 
the darkest curve at the bottom, to W[—0.15,0.15] for the lightest curve at 
the top, in 11 increments. For every jitter level, we plotted the mean relative 
en'or ||D — Z)||ir/|jZi>||ir for all algorithms. 


7, we assume that the observed entries are known exactly, 
and we plot the success rate (percentage of accurate EDM 
reconstructions) against the number of deletions in the first 
case, and the number of calibration events in the second case. 
Accurate reconstruction is defined in terms of the relative error. 
Let D be the true, and ID the estimated EDM. The relative 
error is then \\D — D\\f/\\D\\f, and we declare success if 
this error is below 1%. 

To generate Eigs. 8 and 9 we varied the amount of random, 
uniformly distributed jitter added to the distances, and for each 
jitter level we plotted the relative error. The exact values of 
intermediate curves are less important than the curves for the 
smallest and the largest jitter, and the overall shape of the 
ensemble. 

A number of observations can be made about the perfor¬ 
mance of algorithms. Notably, OptSpace (Algorithm 3) does 
not perform well for randomly deleted entries when n = 20; 
it was designed for larger matrices. Eor this matrix size, 
the mean relative reconstruction error achieved by OptSpace 
is the worst of all algorithms (Eig. 8). In fact, the relative 
error in the noiseless case was rarely below the success 
threshold (set to 1%) so we omitted the corresponding near¬ 
zero curve from Fig. 6. Furthermore, OptSpace assumes that 
the pattern of missing entries is random; in the case of a 
blocked deterministic structure associated with MDU, it never 
yields a satisfactory completion. 

On the other hand, when the unobserved entries are ran¬ 
domly scattered in the matrix, and the matrix is large—in the 
ultrasonic calibration example the number of sensors n was 
200 or more—OptSpace is a very fast and attractive algorithm. 
To fully exploit OptSpace, n should be even larger, in the 
thousands or tens of thousands. 

SDR (Algorithm 5) performs well in all scenarios. For both 


Fig. 9. Comparison of different algorithms applied to multidimensional 
unfolding with varying number of acoustic events k and noisy distances. For 
every number of acoustic events, we generated 1000 realizations of m = 20 
microphone locations uniformly at random in a unit cube. In addition to the 
number of acoustic events, we varied the amount of random jitter added to 
the distances, with the same parameters as in Fig. 8. For every jitter level, we 
plotted the mean relative error ||Z9 — Z9||ir/|jii)||i? for all algorithms. 

the random deletions and the MDU, it has the highest success 
rate, and it behaves well with respect to noise. Alternating 
coordinate descent (Algorithm 4) performs slightly better in 
noise for small number of deletions and large number of 
calibration events, but Figs. 6 and 7 indicate that for certain 
realizations of the point set it gives large errors. If the worst- 
case performance is critical, SDR is a better choice. We note 
that in the experiments involving the SDR, we have set the 
multiplier A in 36 to the square root of the number of missing 
entries. This choice was empirically found to perform well. 

The main drawback of SDR is speed; it is the slowest among 
the tested algorithms. To solve the semidefinite program we 
used CVX [50], [51], a Matlab interface to various interior 
point methods. For larger matrices (e.g., n = 1000), CVX runs 
out of memory on a desktop computer, and essentially never 
finishes. Matlab implementations of alternating coordinate 
descent, rank alternation (Algorithm 2), and OptSpace are all 
much faster. 

The microphone calibration algorithm by Crocco [27] per¬ 
forms equally well for any number of acoustic events. This 
may be explained by the fact that it always reduces the 
problem to ten unknowns. It is an attractive choice for practical 
calibration problems with a smaller number of calibration 
events. Algorithm’s success rate can be further improved if 
one is prepared to run it for many random initializations of 
the non-linear optimization step. 

Interesting behavior can be observed for the rank alternation 
in MDU. Figs. 7 and 9 both show that at low noise levels, the 
performance of the rank alternation becomes worse with the 
number of acoustic events. At first glance, this may seem coun¬ 
terintuitive, as more acoustic events means more information; 
one could simply ignore some of them, and perform at least 
equally well as with fewer events. But this reasoning presumes 
that the method is aware of the geometrical meaning of the 
matrix entries; on the contrary, rank alternation is using only 
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rank. Therefore, even if the percentage of the observed matrix 
entries grows until a certain point, the size of the structured 
blocks of unknown entries grows as well (and the percentage 
of known entries in columns/rows corresponding to acoustic 
events decreases). This makes it harder for a method that does 
not use geometric relationships to complete the matrix. A loose 
comparison can be made to image inpainting: If the pixels are 
missing randomly, many methods will do a good job; but if a 
large patch is missing, we cannot do much without additional 
structure (in our case geometry), no matter how large the rest 
of the image is. 

To summarize, for smaller matrices the SDR seems to be 
the best overall choice. For large matrices the SDR becomes 
too slow and one should turn to alternating coordinate descent, 
rank alternation or OptSpace. Rank alternation is the simplest 
algorithm, but alternating coordinate descent performs better. 
For very large matrices (n on the order of thousands or tens 
of thousands), OptSpace becomes the most attractive solution. 
We note that we deliberately refrained from making detailed 
running time comparisons, due to the diverse implementations 
of the algorithms. 

E Summary 

In this section we discussed: 

• The problem statement for EDM completion and de- 
noising, and how to easily exploit the rank property 
(Algorithm 2), 

• Standard objective functions in MDS: raw stress and 
s-stress, and a simple algorithm to minimize s-stress 
(Algorithm 4), 

• Different semidehnite relaxations that exploit the connec¬ 
tion between EDMs and PSD matrices, 

• Multidimensional unfolding, and how to solve it effi¬ 
ciently using EDM completion, 

• Performance of the introduced algorithms in two very 
different scenarios: EDM completion with randomly un¬ 
observed entries, and EDM completion with deterministic 
block structure of unobserved entries (MDU). 

IV. Unlabeled Distances 

In certain applications we can measure the distances be¬ 
tween the points, but we do not know the correct labeling. 
That is, we know all the entries of an EDM, but we do not 
know how to arrange them in the matrix. As illustrated in Eig. 
10 A, we can imagine having a set of sticks of various lengths. 
The task is to work out the correct way to connect the ends 
of different sticks so that no stick is left hanging open-ended. 

In this section we exploit the fact that in many cases, 
distance labeling is not essential. Eor most point conhgura- 
tions, there is no other set of points that can generate the 
corresponding set of distances, up to a rigid transformation. 

Localization from unlabeled distances is relevant in various 
calibration scenarios where we cannot tell apart distance 
measurements belonging to different points in space. This can 
occur when we measure times of arrivals of echoes, which 
correspond to distances between the microphones and the im¬ 
age sources (see Eig. 12) [29], [6]. Somewhat surprisingly, the 


same problem of unlabeled distances appears in sparse phase 
retrieval; to see how, take a look at the “Phase Retrieval” 
box. 

No efficient algorithm currently exists for localization from 
unlabeled distances in the general case of noisy distances. We 
should mention, however, a recent polynomial-time algorithm 
(albeit of a high degree) by Gujarathi and et al. [31], that can 
reconstruct relatively large point sets from unordered, noiseless 
distance data. 

At any rate, the number of assignments to test is sometimes 
sufficiently small so that an exhaustive search does not present 
a problem. We can then use EDMs to hnd the best labeling. 
The key to the unknown permutation problem is the following 
fact. 

Theorem 3. Draw Xi,X 2 ,--‘ ,£c„ G independently from 
some absolutely continuous probability distribution (e.g. uni¬ 
formly at random) on D, C Then with probability 1, 
the obtained point configuration is the unique (up to a rigid 
transformation) point configuration in Q that generates the set 
of distances {\\xi — Xj\\ ,1 < i < j < n\. 

This fact is a simple consequence of a result by Boutin 
and Kemper [52] who give a characterization of point sets 
reconstructible from unlabeled distances. 

Eigs. 1013 and IOC show two possible arrangements of the 
set of distances in a tentative EDM; the only difference is that 
the two hatched entries are swapped. But this simple swap 
is not harmless: there is no way to attach the last stick in 
Eig. lOD, while keeping the remaining triangles consistent. We 
could do it in a higher embedding dimension, but we insist on 
realizing it in the plane. 

What Theorem 3 does not tell us is how to identify the 
correct labeling. But we know that for most sets of distances, 
only one (correct!) permutation can be realized in the given 
embedding dimension. Of course, if all the labelings are 
unknown and we have no good heuristics to trim the solution 
space, finding the correct labeling is difficult, as noted in [31]. 
Yet there are interesting situations where this search is feasible 
because we can augment the EDM point by point. We describe 
one such situation in the next subsection. 

A. blearing the Shape of a Room [6] 

An important application of EDMs with unlabeled distances 
is the reconstruction of the room shape from echoes. An 
acoustic setup is shown in Eig. 12A, but one could also use 
radio signals. Microphones pick up the convolution of the 
sound emitted by the loudspeaker with the room impulse 
response (RIR), which can be estimated by knowing the 
emitted sound. An example RIR recorded by one of the 
microphones is illustrated in Eig. 12B, with peaks highlighted 
in green. Some of these peaks are hrst-order echoes coming 
from different walls, and some are higher-order echoes or just 
noise. 

Echoes are linked to the room geometry by the image source 
model [53]. According to this model, we can replace echoes by 
image sources (IS)—mirror images of the true sources across 
the corresponding walls. Position of the image source of s 
corresponding to wall i is computed as 
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EDM Perspective on Sparse Phase Retrieval (The Unexpected Distance Structure) 


In many cases, it is easier to measure a signal in the Fourier 
domain. Unfortunately, it is common in these scenarios that we 
can only reliably measure the magnitude of the Fourier transform 
(FT). We would like to recover the signal of interest from just 
the magnitude of its FT, hence the name phase retrieval. X-ray 
crystallography [54] and speckle imaging in astronomy [55] are 
classic examples of phase retrieval problems. In both of these 
applications the signal is spatially sparse. We can model it as 


f{x) = '^CiS{x-Xi), (47) 

where a are the amplitudes and Xi the locations of the n 
Dirac deltas in the signal. In what follows, we discuss the 
problem on 1-dimensional domains, that is a: G R, knowing 
that a multidimensional phase retrieval problem can be solved 
by solving many 1-dimensional problems [7]. 

Note that measuring the magnitude of the FT of f{x) is 
equivalent to measuring its autocorrelation function (ACF). For 
a sparse f{x), the ACF is also sparse and given as 

n n 

a{x) = '^'^CiCjS{x — {xi — Xj)), (48) 

i=i j=i 

where we note the presence of differences between the locations 
Xi in the support of the ACF. As a{x) is symmetric, we do 
not know the order of Xi, and so we can only know these 
differences up to a sign, which is equivalent to knowing the 
distances ||a3i — Xj\\. 

For the following reasons, we focus on the recovery of the support 
of the signal f{x) from the support of the ACF a{x)'. i) in certain 
applications, the amplitudes d may be all equal, thus limiting 


their role in the reconstruction; ii) knowing the support of fix) 
and its ACF is sufficient to exactly recover the signal f{x) [7]. 
The recovery of the support of fix) from the one of a(cc) 
corresponds to the localization of a set of n points from their 
unlabeled distances: we have access to all the pairwise distances 
but we do not know which pair of points corresponds to any 
given distance. This can be recognized as an instance of the 
turnpike problem, whose computational complexity is believed 
not to be NP-hard but for which no polynomial time algorithm 
is known [56]. 

From an EDM perspective, we can design a reconstruction 
algorithm recovering the support of the signal fix) by labeling 
the distances obtained from the ACF such that the resulting 
EDM has rank smaller or equal than 3. This can be regarded 
as unidimensional scaling with unlabeled distances, and the 
algorithm to solve it is similar to echo sorting (Algorithm 6). 



Fig. 11. A graphical representation of the phase retrieval problem for 
1-dimensional sparse signals. (A) We measure the ACF of the signal 
and we recover a set of distances (st/cks in Fig. 10) from its support. 
(B) These are the unlabeled distances between all the pairs of Dirac 
deltas in the signal fix). We exactly recover the support of the signal 
if we correctly label the distances. 


Si = s -F 2 (p.j - s, rii) rii, (49) 

where Pi is any point on the ith wall, and rii is the unit 
normal vector associated with the fth wall, see Fig. 12A. A 
convex room with planar walls is completely determined by 
the locations of first-order ISs [6], so by reconstructing their 
locations, we actually reconstruct the room’s geometry. 

We assume that the loudspeaker and the microphones are 
synchronized so that the times at which the echoes arrive 
directly correspond to distances. The challenge is that the 
distances—the green peaks in Fig. 12B— are unlabeled; it 
might happen that the kth peak in the RIR from microphone 
1 and the A:th peak in the RIR from microphone 2 come from 
different walls, especially for larger microphone arrays. Thus, 
we have to address the problem of echo sorting, in order to 
group peaks corresponding to the same image source in RIRs 
from different microphones. 

Assuming that we know the pairwise distances between the 
microphones R— [ri,..., r^], we can create an EDM corre¬ 
sponding to the microphone array. Because echoes correspond 
to image sources, and image sources are just points in space, 
we attempt to grow that EDM by adding one point—an image 
source—at a time. To do that, we pick one echo from every 
microphone’s impulse response, augment the EDM based on 
echo arrival times, and check how far the augmented matrix is 


from an EDM with embedding dimension three, as we work in 
3D space. The distance from an EDM is measured with s-stress 
cost function. It was shown in [6] that a variant of Theorem 
3 applies to image sources when microphones are thrown at 
random. Therefore, if the augmented matrix satisfies the EDM 
properties, almost surely we have found a good image source. 
With probability 1, no other combination of points could have 
generated the used distances. 

The main reason for using EDMs and s-stress instead of, for 
instance, the rank property, is that we get robust algorithms. 
Echo arrival times are corrupted with various errors, and 
relying on the rank is too brittle. It was verified experimentally 
[6] that EDMs and s-stress yield a very robust filter for the 
correct combinations of echoes. 

Thus we may try all feasible combinations of echoes, and 
expect to get exactly one “good” combination for every image 
source that is “visible” in the impulse responses. In this case, 
as we are only adding a single point, the search space is 
small enough to be rapidly traversed exhaustively. Geometric 
considerations allow for a further trimming of the search 
space: because we know the diameter of the microphone array, 
we know that an echo from a particular wall must arrive at all 
the microphones within a temporal window corresponding to 
the array’s diameter. 

The procedure is as follows; collect all echo arrival times 
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Fig. 10. Illustration of the uniqueness of EDMs for unlabeled distances. A 
set of unlabeled distance (A) is distributed in two different ways in a tentative 
EDM with embedding dimension 2 (B and C). The correct assignment yields 
the matrix with the expected rank (C), and the point set is easily realized in 
the plane (E). On the contrary, swapping just two distances (hatched squares 
in (B) and (C)) makes it impossible to realize the point set in the plane (D). 
Triangles that do not coincide with the swapped edges can still be placed, but 
in the end we are left with a hanging orange stick that cannot attach itself to 
any of the five nodes. 


received by the jth microphone in the set Ti, and fix ti S Ti 
corresponding to a particular image source. Then Algorithm 
6 finds echoes in other microphones’ RlRs that correspond 
to this same image source. Once we group all the peaks 
corresponding to one image source, we can determine its 
location by multilateration {e.g. by running the classical MDS), 
and then repeat the process for other echoes in Ti. 

To get a ballpark idea of the number of combinations to test, 
suppose that we detect 20 echoes per microphone^, and that the 
diameter of the five-microphone array is 1.5 m. Thus for every 
peak time ti G Ti we have to look for peaks in the remaining 
four microphones that arrived within a window around ti of 
length 2 x where 343 m/s is the speed of sound. This 

is approximately 6 ms, and in a typical room we can expect 
about five early echoes within a window of that duration. Thus 
we have to compute the s-stress for 20 x 5“^ = 12500 matrices 
of size 6x6, which can be done in a matter of seconds on a 
desktop computer. In fact, once we assign an echo to an image 
source, we can exclude it from further testing, so the number 
of combinations can be further reduced. 

Algorithm 6 was used to reconstruct rooms with centimeter 
precision [6] with one loudspeaker and an array of five micro¬ 
phones. The same algorithm also enables a dual application: 
indoor localization of an acoustic source using only one 


^We do not need to look beyond early echoes corresponding to at most 
three bounces. This is convenient, as echoes of higher orders are challenging 
or impossible to isolate. 


Fig. 12. (A) Illustration of the image source model for first- and second- 

order echoes. Vector rii is the outward-pointing unit normal associated with 
the ith wall. Stars denote the image sources, and Sij is the image source 
corresponding to the second-order echo. Sound rays corresponding to first 
reflections are shown in blue, and the ray corresponding to the second-order 
reflection is shown in red. (B) Early part of a typical recorded room impulse 
response. 


Algorithm 6 Echo Sorting [6] 


function ECHOSORT(i2, tl, ■ ■ ■ , Tm) 

D G- edm(i2) 

Sbest -flnf 

for all t = [t2, ■ ■ ■, tm], such that ti G Ti do 

d •(— c • [tl, [> c is the sound speed 

d\ 

dJ 0 


Z?aug ^ 


7: if s—stress (Daug) < Sbest then 

8: Sbest ^ s-stress (Daug) 

9: dbest G- d 

10 : end if 

11 : end for 

12 : return dbest 

13: end function 


microphone—a feat not possible if we are not in a room [57]. 


B. Summary 

To summarize this section: 

• We explained that for most point sets, the distances 
they generate are unique; there are no other point sets 
generating the same distances, 

• In room reconstruction from echoes, we need to identify 
the correct labeling of the distances to image sources. 
EDMs act as a robust filter for echoes coming from the 
same image source. 
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TABLE II 

Applications of EDMs with different twists 


Application 

missing 

distances 

noisy 

distances 

unlabeled 

distances 

Wireless sensor networks 

✓ 

✓ 

X 

Molecular conformation 

✓ 

✓ 

X 

Hearing the shape of a room 

X 

✓ 

✓ 

Indoor localization 

X 

✓ 

✓ 

Calibration 

✓ 

✓ 

X 

Sparse phase retrieval 

X 

✓ 

✓ 


• Sparse phase retrieval can be cast as a distance problem, 
too. The support of the ACF gives us distances between 
the deltas in the original signal. Echo sorting can be 
adapted to solve the problem from the EDM perspective. 

V. Ideas for Euture Research 

Even problems that at first glance seem to have little to do 
with EDMs, sometimes reveal a distance structure when you 
look closely. A good example is sparse phase retrieval. 

The purpose of this paper is to convince the reader that 
Euclidean distance matrices are powerful objects with a mul¬ 
titude of applications (Table 11 lists various flavors), and that 
they should belong to any practitioner’s toolbox. We have 
an impression that the power of EDMs and the associated 
algorithms has not been sufficiently recognized in the signal 
processing community, and our goal is to provide a good 
starting reference. To this end, and perhaps to inspire new 
research directions, we list several EDM-related problems that 
we are curious about and believe are important. 

a) Distance matrices on manifolds: If the points lie on 
a particular manifold, what can be said about their distance 
matrix? We know that if the points are on a circle, the 
EDM has rank three instead of four, and this generalizes to 
hyperspheres [17]. But what about more general manifolds? 
Are there invertible transforms of the data or of the Gram 
matrix that yield EDMs with a lower rank than the embedding 
dimension suggests? What about different distances, e.g. the 
geodesic distance on the manifold? Answers to these questions 
have immediate applications in machine learning, where the 
data can be approximately assumed to be on a smooth surface 
[23]. 

b) Projections of EDMs on lower dimensional subspaces: 
What happens to an EDM when we project its generating 
points to a lower dimensional space? What is the minimum 
number of projections that we need to be able to recon¬ 
struct the original point set? Answers to these questions have 
significant impact on imaging applications such as X-ray 
crystallography and seismic imaging. What happens when we 
only have partial distance observations in various subspaces? 
What are other useful low-dimensional structures on which we 
can observe the high-dimensional distance data? 

c) Efficient algorithms for distance labeling: Without 
application-specific heuristics to trim down the search space, 
identifying correct labeling of the distances quickly becomes 
an arduous task. Can we identify scenarios for which there 
are efficient labeling algorithms? What happens when we do 


not have the labeling, but we also do not have the complete 
collection of sticks! What can we say about uniqueness of 
incomplete unlabeled distance sets? Some of the questions 
have been answered by Gujarathi [31], but many remain. The 
quest is on for faster algorithms, as well as algorithms that 
can handle noisy distances. 

In particular, if the noise distribution on the unlabeled 
distances is known, what can we say about the distribution 
of the reconstructed point set (taking in some sense the best 
reconstruction over all labelings)? Is it compact, or we can 
Jump to totally wrong assignments with positive probability? 

d) Analytical local minimum of s-stress: Everyone agrees 
that there are many, but to the best of our knowledge, no 
analytical minimum of s-stress has yet been found. 

VI. Conclusion 

At the end of this tutorial, we hope that we succeeded 
in showing how universally useful EDMs are, and that we 
inspired readers coming across this material for the first time 
to dig deeper. Distance measurements are so common that a 
simple, yet sophisticated tool like EDMs deserves attention. 
A good example is the semidefinite relaxation; even though it 
is generic, it is the best performing algorithm for the specific 
problem of ad-hoc microphone array localization. Continuing 
research on this topic will bring new revolutions, like it did in 
the 80s in crystallography. Perhaps the next one will be fueled 
by solving the labeling problem. 
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